Part 1: Spatial autocorrelation in GeoDa

Question 1: what is the Global Moran’s I and what does it suggest about the variable? Repeat for the other two weight matrices.

Weights matrix used Global Moran’s I value
contiguity, 1st order 0.472
contiguity, 2nd order 0.212
distance-based 0.283

All of these Moran’s I values are postive, indicating that the variable shows positive spatial autocorrelation. The Moran’s I value changes depending on the weights matrix used. The Moran’s I value calculated using the 1st order contiguity weights matrix is the most postive. This weights matrix considers all contiguous polygons using queens adjacency to be neighbors. The 2nd order contiguity takes neighbors of neighbors, so it makes sense that using this weights matrix we see less strong positive autocorrelation. The distance based weights matrix takes the distance between centroids to determine adjacency – this ends up giving fewer neighbors to larger polygons. This may have something to do with its relatively small Moran’s I value compared to the contiguity, 1st order value.

Question 2a: for each weight matrix, how many districts have Moran’s I values that are significant at p=0.001, and how many districts are in high-high (hot-spot) clusters and low-low (clod-spot) clusters.

Weights matrix used # significant Moran’s I # in high-high clusters # in low-low clusters
contiguity, 1st order 7 16 8
contiguity, 2nd order 9 17 10
distance-based 8 11 13

Question 2b: Repeat for local G cluster map

Weights matrix used # significant Gi # in high clusters # in low clusters
contiguity, 1st order 7 17 9
contiguity, 2nd order 7 19 16
distance-based 7 11 17

Spatial auto-correlation in R

Question 1. create a histogram of the values for MALKIDS [hint: hist()]

library(rgdal)
## Loading required package: sp
## rgdal: version: 1.2-5, (SVN revision 648)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-4
datapath = "/users/miadawson/Documents/~UC Davis~/~GEO 200CN/labs/lab 9"
nepal=readOGR(dsn=datapath, layer="Nepal")
## OGR data source with driver: ESRI Shapefile 
## Source: "/users/miadawson/Documents/~UC Davis~/~GEO 200CN/labs/lab 9", layer: "Nepal"
## with 75 features
## It has 26 fields
## Integer64 fields read as strings:  ID KIDS1_5 AD_ILGT50
hist(nepal$MALKIDS, main='Rate of childhood malnutrition in Nepalese districts', xlab = NULL)

Question 2: let’s also map the values, but now that we’ve been spoiled by the interactive maps in GeoDa let’s try to do something similar using the tools in the ggplot2 and tmap packages, designed to make more complex graphs and maps than the base graphing package.

library (tmap)
#vignette(package = "tmap") # available vignettes in tmap
#vignette("tmap-nutshell") #gives you an overview of the package
tmap_mode('plot')
## tmap mode set to plotting
qtm (shp=nepal, fill="MALKIDS", fill.palette="Blues") #map just the MALKIDS

qtm (shp=nepal, fill=c("MALKIDS", "POPULATION"), fill.palette="Blues") #map both

tmap has two modes “plot” return static maps and “view” returns interactive maps. You toggle between them with the function tmap_mode() [ tmap_mode(”view”) or tmap_mode(”plot”)]. Switch to “view” mode and plot the maps above again.

tmap_mode('view')
## tmap mode set to interactive viewing
qtm (shp=nepal, fill="MALKIDS", fill.palette="Blues") #map just the MALKIDS
qtm (shp=nepal, fill=c("MALKIDS", "POPULATION"), fill.palette="Blues") #map both

Now let’s plot a scatterplot similar to what we saw in GeoDa.

library (ggplot2)
library(labeling)
ggplot(nepal@data, aes(POPULATION, MALKIDS)) + geom_point(aes(col=MALKIDS, size=MALKIDS))

Both tmap and ggplot2 work by adding layers so in this case we set the general plot tell it that the data comes from npl@data, and the values for the aesthetic mapping come from population, and malkids. Then we add another layer, the actual points where the color and size of the points are defined by the value of malkids.

Looking at the polt you might decide you want the higher values to have the darker colors, you can do that by simply adding a negative sign in front of MALKIDS [col=–MALKIDS].

Review the material from Lab 8 to:

1) create a) 1st order queen contiguity; b) 2nd order queen contiguity; and c) distance based (0,0,7) weight matrices.

library(raster)
library(spdep)
## Loading required package: Matrix
##1st order queen cont
?poly2nb
w <- poly2nb(nepal, row.names=nepal$ID) ##default is queen=TRUE
ww <-  nb2listw(w, style='B')

##2nd order queen cont
w2 = nblag(w, 2)
w22=w2[[2]] #subset for only 2nd order neighbors
ww22 <-  nb2listw(w22, style='B')

##distance-based cont
dw = poly2nb(nepal, row.names=nepal$ID, snap = 0.7) 
dww <-  nb2listw(dw, style='B')

2) calculate the global Moran’s I with each weight matrix.

##1st order queen cont
moran(nepal$MALKIDS, ww, n=length(ww$neighbours), S0=Szero(ww))
## $I
## [1] 0.4798273
## 
## $K
## [1] 2.430716
##2nd order queen cont
moran(nepal$MALKIDS, ww22, n=length(ww22$neighbours), S0=Szero(ww22)) 
## $I
## [1] 0.2258045
## 
## $K
## [1] 2.430716
##distance-based cont
moran(nepal$MALKIDS, dww, n=length(dww$neighbours), S0=Szero(dww)) 
## $I
## [1] 0.2616114
## 
## $K
## [1] 2.430716

Are the values the same as those you got in GeoDa, what might drive variation?

Weights matrix used Moran’s I value (GeoDa) Moran’s I Value (R)
contiguity, 1st order 0.472 0.479
contiguity, 2nd order 0.212 0.225
distance-based 0.283 0.261

These values are close, but not exactly the same as the values I got in GeoDa. There may be slight differences in the way GeoDa and the R package calculate these values. For example, if the programs round differently within a complicated expression, this could lead to different output values for Moran’s I.

Now let’s move into looking at some of the local patterns. First we’ll calculate the Local Moran’s I, and create a dataframe with the values, lag values and p-vlaues to then we visulaize the results in similar ways to what we saw in GeoDa.

Again test with each of your spatial weight matrcies.

##1st order queens
LoM <- localmoran(nepal$MALKIDS, ww, zero.policy=NULL)
LoMdf <- data.frame (LoM)
lag = lag.listw(ww, nepal$MALKIDS)
dat = data.frame(ID = nepal$ID, mal = nepal$MALKIDS, lag = lag, pval = LoMdf$Pr.z...0)
dat$sig =ifelse (dat$pval < 0.05 ,1, 0 ) #set signifcance at p=0.05

mv = mean(dat$mal) #create a mean for the value
mh = mean (dat$lag) #create a mean for the lagged value

ggplot(dat, aes(x = mal, y = lag, col=-sig)) + geom_point() + geom_smooth(method = "lm", col="grey") +
geom_vline(xintercept=mv, linetype = "dashed", col="grey") + geom_hline(yintercept = mh, linetype="dashed",
col="grey")

##2nd order queens
LoM <- localmoran(nepal$MALKIDS, ww22, zero.policy=NULL)
LoMdf <- data.frame (LoM)
lag = lag.listw(ww22, nepal$MALKIDS)
dat = data.frame(ID = nepal$ID, mal = nepal$MALKIDS, lag = lag, pval = LoMdf$Pr.z...0)
dat$sig =ifelse (dat$pval < 0.05 ,1, 0 ) #set signifcance at p=0.05

mv = mean(dat$mal) #create a mean for the value
mh = mean (dat$lag) #create a mean for the lagged value

ggplot(dat, aes(x = mal, y = lag, col=-sig)) + geom_point() + geom_smooth(method = "lm", col="grey") +
geom_vline(xintercept=mv, linetype = "dashed", col="grey") + geom_hline(yintercept = mh, linetype="dashed",
col="grey")

##distance-based
LoM <- localmoran(nepal$MALKIDS, dww, zero.policy=NULL)
LoMdf <- data.frame (LoM)
lag = lag.listw(dww, nepal$MALKIDS)
dat = data.frame(ID = nepal$ID, mal = nepal$MALKIDS, lag = lag, pval = LoMdf$Pr.z...0)
dat$sig =ifelse (dat$pval < 0.05 ,1, 0 ) #set signifcance at p=0.05

mv = mean(dat$mal) #create a mean for the value
mh = mean (dat$lag) #create a mean for the lagged value

ggplot(dat, aes(x = mal, y = lag, col=-sig)) + #sets x axis as malnutrition rate, y value as lagged value, color as significance 
  geom_point() + #specifies scatterplot
  geom_smooth(method = "lm", col="grey") + #adds smoothed line of best fit
  geom_vline(xintercept=mv, linetype = "dashed", col="grey") + #adds vertical line bisecting data
  geom_hline(yintercept = mh, linetype="dashed", col="grey") #adds horizontal line bisecting data

Try to unpick this rather long line, and see if you understand the different layers that are making up this plot. Looking at this plot what can you say about the pattern of Local Moran’s I? (see above comments)

Looking at this plot what can you say about the pattern of Local Moran’s I?

From these plots, it appears that there is most postive spatial autocorrelation when using the 1st order queens contiguity. According to OSU, positive spatial autocorrelation is evidenced by more data falling in the bottom left and upper right corners of Moran’s scatterplots. This pattern is most visibile in the plot using 1st order quuen’s contiguity. So, it should have the highest local Moran’s I. The 2nd order queen’s contiguity based plot appears to have a less positive local moran’s I, and the distance based seems to have the least positive local moran’s I.

For the Gi and Gi* go read and follow along with the with the Rspatial chapter on local Statistics to create local statistics for the Nepal dataset[http://rspatial.org/rosu/rst/Chapter8.html], stop with the Local Moran’s I ie don’t work through the Geographically Weighted Regression section.

##local getis Gi
wr <- poly2nb(nepal, row.names=nepal$ID, queen=FALSE)
lstw <- nb2listw(wr, style='B')
Gi <- localG(nepal$MALKIDS, lstw)
head(Gi)
## [1]  0.4309344  1.5261740  0.2739212 -0.8286501 -3.2037288 -1.7890001
#plot 
par(mai=c(0,0,0,0))
Gcuts <- cut(Gi, 5)
Gcutsi <- as.integer(Gcuts)
cols <- rev(gray(seq(0,1,.2)))
plot(nepal, col=cols[Gcutsi])
legend('bottomleft', levels(Gcuts), fill=cols)

##classify each polygon as own neighbor
ws <- include.self(wr)
lstws <- nb2listw(ws, style='B')
Gis <- localG(nepal$MALKIDS, lstws)
Gscuts <- cut(Gis, 5)
Gscutsi <- as.integer(Gscuts)
cols <- rev(gray(seq(0,1,.2)))
plot(nepal, col=cols[Gscutsi])
legend('bottomleft', levels(Gscuts), fill=cols)

#local average
m <- sapply(ws, function(i) mean(nepal$MALKIDS[i]))
cts <- cut(m, 5)
mcts <- as.integer(cts)
plot(nepal, col=cols[mcts])
legend('bottomleft', levels(cts), fill=cols)

#local Moran's i
Ii <- localmoran(nepal$MALKIDS, lstw)
Icuts <- cut(Ii, 5)
Icutsi <- as.integer(Icuts)
plot(nepal, col=cols[Icutsi])
legend('bottomleft', levels(Icuts), fill=cols)

Again, where the values different from those you got with GeoDa, and what might drive variation?

These values are similar to what I got in GeoDa. As before, there may be differences in the way GeoDa and the R package calculate these values. For example, if the programs round differently within a complicated expression, this could lead to different output values.